clear
rng default
addpath /Users/cristinagualdani/mosek/9.2/toolbox/r2015aom
%% Load useful matrices
load PYX_cond 
load U_CS

%% Construct grid of parameters to check for each realisation of X
u0grid_temp=0; %[Ychosen=0, X=1] [Ychosen=0, X=2], location normalisation 
u1grid_temp_all=[U_CS(1,1) U_CS(2,1) U_CS(3,1)]; %[Ychosen=1, X=1] [Ychosen=1, X=2] [Ychosen=1, X=3], scale normalisation
gran=0.1; 
u2grid_temp=-20:gran:20;  %[Ychosen=2, X=1] [Ychosen=2, X=2] [Ychosen=2, X=3]
u3grid_temp=-20:gran:20;  %[Ychosen=3, X=1] [Ychosen=3, X=2] [Ychosen=3, X=3]
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)

%% Construct identified set for each realisation of X
ntypes_m=3;
IdSetM=cell(ntypes_m,1);


for x=1:ntypes_m %we can proceed separately across types as long as we do not impose independence from X
    u1grid_temp=u1grid_temp_all(x);
    [ca, cb, cc,cd] = ndgrid(u0grid_temp, u1grid_temp, u2grid_temp, u3grid_temp);
    u0grid=ca(:);
    u1grid=cb(:);
    u2grid=cc(:);
    u3grid=cd(:);
    sg=size(u0grid,1);
    % Components invariant across parameter values
    [sgreduced, ...
    id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
    P_col,...
    Tcoord_zeros,...
    Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
    cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros,...
    equalorder, n_vertices,infty_pos]=useful_anyparam(u0grid, u1grid, u2grid, u3grid, sg);
    % Loop over parameter values, order of columns: [u0x u1x u2x u3x]
    IdSetM{x}=OneSide(u0grid, u1grid, u2grid, u3grid, ...
                       sgreduced, PYX_cond(x,:), ...
                       id_Ugridreduced, s_id_gridreduced, sU_id_gridreduced,...
                       P_col, ...
                       Tcoord_zeros,...
                       Ugridreduced, n_U_sets, indices_pairs_extended, Tcoord, C, C_sign,...
                       cr_obsequiv, cr_one, cr_zeros,fill_obsequiv, fill_one, fill_zeros, numeqconstraints, RHS_one, RHS_zeros, equalorder,n_vertices,infty_pos);  
end

save('IdSetM.mat', 'IdSetM')

